Self-induced density modulations in the free expansion of Bose-Einstein condensates 
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We simulate numerically the free expansion of a repulsive Bose-Einstein condensate with an ini- 
tially Gaussian density profile. We find a self-similar expansion only for weak inter-atomic repulsion. 
In contrast, for strong repulsion we observe the spontaneous formation of a shock wave at the sur- 
face followed by a significant depletion inside the cloud. In the expansion, contrary to the case of 
a classical viscous gas, the quantum fluid can generate radial rarefaction density waves with sev- 
eral minima and maxima. These intriguing nonlinear effects, never observed yet in free-expansion 
experiments with ultra-cold alkali-metal atoms, can be detected with the available setups. 
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The anisotropic free expansion of a gas of 87 Rb atoms 
was the first experimental evidence of Bose-Einstein con- 
densation in ultra-cold gases The non-ballistic free 
expansion observed with 6 Li atoms has been saluted 
as the first signature of superfiuid behavior in a ultra- 
cold fermi vapor Q. In both cases the atomic quantum 
gases can be described by the hydrodynamic equations 
of superfiuids and, because the initial density profile is 
an inverted parabola, the free expansion is self-similar 

The first theoretical investigations of the free expan- 
sion into vacuum of a classical gas sphere with constant 
initial density dates back to the 1960's numerical 
■ analysis was needed to analyze in detail the formation 
of a depletion at the center and of a shock wave at the 
surface [ToL 11 1. More recently, rarefaction waves have 
been produced experimentally by the free expansion of 
an electron plasma 12 1 . 

In this Letter we show that the free expansion of a 
bosonic superfiuid into vacuum displays intriguing non- 
linear phenomena. In particular, by integrating numer- 
ically [13| the time-dependent Gross-Pitaevskii equation 
[f"H we prove that the expansion of a repulsive Bose- 
Einstein condensate (BEC) of initial Gaussian density 
profile gives rise to self-induced density modulations, i.e. 
self-induced rarefaction waves. In addition, we find that 
the bosonic cloud produces a shock wave at the surface, 
that is damped by the quantum pressure of the super- 
fluid. These nonlinear effects, which can be observed ex- 
perimentally with available techniques, are strongly sup- 
pressed if a harmonic confinement along two perpendic- 
ular directions is retained and free expansion is allowed 
in ID only. Observe that the formation of shock and rar- 
efaction waves induced by an ad hoc external perturba- 



tion has been suggested 15, 1(| 17 1 and observed recently 
m BECs [3 [3 

and also n non-dissipative nonlinear op- 
tics 

We describe the collective motion of the BEC of N 
atoms in terms of a complex mean-field wavefunction 
i]j(r,t) normalized to unity and such that p(r,t) — 
Npi(r,t) — N \yj(r, t)\ 2 is the number-density distribu- 
tion. The equation of motion that we assume for ip(r,t) 
is the time-dependent Gross-Pitaevskii equation (GPE) 
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where U(r,t) is a confining potential that we assume to 
vanish at t > (thus allowing free expansion) and m 
is the atomic mass. The nonlinear term represents the 
inter-atomic interaction at a mean-field level, where a s is 
the s-wave scattering length, and we consider the repul- 
sive regime a s > 0. 

In traditional experiments with ultra-cold alkali-metal 
atoms [l|, [3| expansion starts from an initial state co- 
inciding with the ground state of the confined superfiuid 
under the action of a (often anisotropic) harmonic poten- 
tial. For robust interparticle interaction (large number of 
particles) , the density profile in this intial state resembles 
closely a negative-curvature parabola When such 

density profile is taken as the initial state of a succes- 
sive free non-ballistic expansion, to a very good degree of 
approximation it expands in a self-similar fashion, main- 
taining the same shape and only spreading out and scal- 
ing down its height proportionally, until a ballistic regime 
is reached when dilution leads to a fully non-interacting 
regime [1I1SS0. 

In the present work we discuss the much more exciting 
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FIG. 1: (Color on line). Four successive frames of the radial 
density profile (solid line) for the expansion of a strongly inter- 
acting condensate, characterized by dimensionless interaction 
parameter g = 2000. The "opacity" of the expanding cloud 
(dashed line) given by the density pi (r) integrated along lines 
at a distance r from the center (rescaled by a factor 0.036). 
All quantities are dimensionless: lengths in units of the initial 
Gaussian width a and time in units of ma 2 /%. 



phenomena observed in the expansion of an interacting 
condensate starting off as a stationary Gaussian 
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FIG. 2: (Color on line). Characteristic times of the free 
expansion as a function of the interaction strength g. For 
g > 48.7, a first minimum in the radial density profile ap- 
pears at the center r = at time £ am (dashed) and fills in 
at tdm (dot-dashed). For g > 363, a local maximum re-forms 
at the center at time t&M (solid) and disappears at time tdM- 
A density maximum appears within the rarefaction region at 
time i C M (dot-dot-dashed). Inset: the time of appearance ta. 
and disappearance td of a denser ring (sketched in the small 
square for g — 1200, t = 1.8) of visibility h (defined in Fig. [TJ 
1% (solid) and 5% (dashed) in the opacity profile (dashed 
lines of Fig. [TJ, as a function of g. Units as in Fig. [T] 



with er 3 = <j x a y <j z . Note that a Gaussian profile is readily 
achieved experimentally by equilibrating the BEC with a 
very small scattering length a s obtained by means of the 



Fano-Feshbach resonace technique [1J, [21( with a care- 
fully tuned external constant magnetic field. a s can then 
be set to the desired value by a sudden change in the 
magnetic field at the time when the harmonic trapping 
potential is switched off. 

The isotropic case a x = a y = a z = a is conceptually 
advantageous, as the expansion can be studied within the 
GPE model in its full generality as a function of a single 
parameter. Consider rescaling the variables of Eq. {T]) as 
follows: r — > r/er, t — * th/(ma 2 ), and ip — » ipa z / 2 ^ to 
produce a dimensionless form 
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of the equation for the free expansion of a conden- 
sate expansion starting off from a Gaussian initial state 
of unit width. The dimensionless interaction strength 
g = 4tt N a s /(T is the one free parameter determining the 
properties of the free expansion. A gas of non interacting 
(g = 0) bosons starting from the Gaussian state expands 



self similarly according to the formula [22j 
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By using an efficient finite-difference Crank-Nicolson al- 
gorithm [l3| we have verified that little changes affect 
the expansion as long as the interaction is small, g < 1 
tracks the weak-coupling limit where interaction only ac- 
celerates slightly the free expansion, while j > 1 repre- 
sents the strong-coupling regime, where the mean-field 
self-interaction term in Eq. §3§ dominates the expansion 
for long enough to produce substantial nonlinear effects 
such as those sketched in Fig. [T] In particular we ob- 
serve the rapid build-up of a sharp expanding spherical 
density wave which leaves behind a central region of de- 
pleted density. New successively formed radial ripples 
cross this density-depleted region. Eventually, at very 
long times, when the overall density has decayed enough 
for the nonlinear term in Eq. ([3]) to become negligible 
everywhere, the expansion recovers a bell-shaped profile. 

More quantitatively, for increasing g, the rarefaction 
starts to be evidenced by a local minimum at the droplet 
center for g > 48.3, and becomes more and more pro- 
nounced and long lived for larger interaction strength g. 
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useful to introduce its local phase velocity, given by 
V= 2 W ' 



(5) 



FIG. 3: (Color on line). Comparison of the density profile 
pi{r) — |?/)(r)| 2 and of the velocity profile v(r)/r for the 
spherically symmetrical expansion of the strongly-interacting 
bosonic cloud. Units as in Fig. [T] 



In the highly nonlinear regime, the free expansion of the 
bosonic cloud into vacuum develops sequences of radial 
density waves with minima and maxima, each starting 
at a characteristic time and disappearing at a later time. 
Figure fj] tracks a few early times in this class, as a func- 
tion of g: more could be defined for larger g. Overall, the 
time i am of appearance of the first local minimum reduces 
slowly as g increases, while the time t^m of disappearance 
of all local minima and recovery of a bell-shaped profile 
increases rapidly with g. In between these two character- 
istic times several traveling density minima and maxima 
can be formed depending on the value of g. Generally, 
the number and density difference of local minima and 
maxima increases with the interaction strength. 

In practice, the visibility of the central depletion need 
not be easy to appreciate by means of total opacity mea- 
surements, since such measurements address the inte- 
grated density of a generic linear section crossing the 
droplet at a given distance from its center. However, 
Fig. [T] shows that for suitably strong interaction, a sen- 
sibly higher-opacity outer ring does indeed develop. The 
inset of Fig. [5] tracks the time range when this inner 
optical-density minimum remains lower than the outside 
denser ring by at least 1% and 5%: it is seen that for 
strong enough interaction (g > 215 and g > 310 respec- 
tively), the depleted region realizes a limited but signifi- 
cant visibility, which improves for stronger coupling. An 
example of visible opacity ring is displayed in the upper 
right corner of Fig. [2] 

To characterize the dynamics of the bosonic cloud it is 



This velocity can be written as v(r, t) = V6>(r, t), where 
8(r, t) is the phase of the macroscopic wave function 
ip(r,t) — p\ir, i) 1 / 2 e l9 ( r '*) of the bosonic supcrfluid. 
From Eq. ([4J one finds immediately the radial phase ve- 
locity of a spherical non-interacting Bosc gas (g = 0): 
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In the interacting case (g > 0) the nonlinear term acts 
as the chemical potential of a fluid of classical pressure 
P = g\ip\ 4 /2 and sound velocity c s = \J g\ip\ 2 [141 ]. For 
g <C 1 the gas velocity follows Eq. ([6|) closely, while 
for large g deviations from Eq. ([6]) are substantial, as 
mainly the interaction term, rather than the quantum 
tendency to derealization, drives the expansion. The ra- 
tio v/r, a constant as a function of r according to Eq. ([6|) 
in a non-interacting expansion, shows strong deviations 
induced by the nonlinear term glV^ 2 - In practice, the 
phase velocity of a strongly interacting BEC approaches 
the local sound velocity, with larger densities implying 
higher velocities. Accordingly, initially the central part 
of the strongly repulsive bosonic superfluid accelerates 
and propagates faster than the periferic part: the en- 
suing mass flow is responsible for the formation of the 
rarefaction region inside the cloud, shown in Fig. [TJ The 
matter flowing quickly out of the central region accumu- 
lates near the profile edge on top of the slower external 
tail, thus tending to produce a shock wave H, 3, 2fjj ]. 
with a BEC density profile extremely steep at the sur- 
face, approaching a step function: this is illustrated by 
the i = 0.24 panel of Fig. [T] This steep wave front survives 
for a brief period, after which density oscillations shoot 
backwards and the surface profile rapidly smoothens its 
density gradient. At this point, the expansion dynam- 
ics is strongly affected by these backward density oscil- 
lations, which induce a inverted relation between local 
velocity and density, as Fig. [3] illustrates: for a strongly 
interacting BEC, at a fixed time t the ratio v/r finds local 
maxima (minima) in correspondence to the local minima 
(maxima) of the density profile p\ . This inversion demon- 
strates the inward motion of the density ripples. These 
local minima represent the rarefaction waves produced 
by the surface step smoothing. This smoothing is driven 
by the quantum pressure term [l6| . 

The quantum pressure — (|V"0| 2 )/(2pi), which plays a 
negligible role in the self-similar non-ballistic expansion 
0] of both Fermi and Bose superfluids with an inverted- 
parabola initial profile [23| , becomes relevant in regular- 
izing the shock-wave singularity, like the dissipative term 
in classical hydrodynamcs . Analogous depletion and 
shock-wave phenomena are indeed observed in the hy- 
drodynamical expansion of hot classical fluids, and can 
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FIG. 4: (Color on line). Self-depleting radial density pro- 
file pi = p/N during the expansion: comparison between the 
strongly interacting Bose gas (GPE, solid lines) and the clas- 
sical gas (NSE, dashed lines), with g — 10 4 and the same 
initial conditions (a unit Gaussian). In the NSE, the shear 
viscosity coefficient r\ = 10~ J . Units as in Fig. [T] 



be simulated e.g. by means of the Navier-Stokes equa- 
tions (NSE), which depend on the (dissipative) coeffi- 
cient of shear viscosity r\ Q. For r/ = the irrotational 
(VAv = 0) NSE reduce to the Euler equations of an ideal 
(non- viscous) fluid. By using P(p) = gPi/2 as the equa- 
tion of state, the Euler equations are exactly equivalent 
to the GPE without the quantum pressure term [3| ■ In 
Fig. g] we compare the expanding BEC (GPE, Eq. fl3J) 
and classical gas (NSE), for g = 10000. The NSE are 
solved by means of a Lagrangian finite-difference method 
[24] ]. The four successive frames displayed in Fig. [4] show 
that, while expanding, the interacting BEC and the clas- 
sical gas produce a remarkably similar self-depletion of 
the central region. On the other hand, the multi-peak 
rarefaction structures predicted by the GPE are not re- 
produced by the NSE. Thus the novelty of the free ex- 
pansion of a BEC with respect to that of a weakly vis- 
cous classical fluid, stands mainly in the distinct large- 
amplitude rarefaction waves moving inside the depletion 
region, which are supported by the nondissipative nature 
of quantum pressure. 

The expansion of the BEC in an anisotropic context 
is also interesting. In particular, we find that a ID 
expansion starting from a Gaussian state produces no 
shock wave and no central depletion. This applies both 
for an idealized ID BEC, represented by a purely ID 
GPE, and for the 3D expansion of a Gaussian wavepacket 
to which an harmonic confining potential is kept along 
two orthogonal space directions. Interestingly, if for the 
initial state of the 3D geometry a spherical Gaussian 
([2]) of width equal to the confining-potcntial harmonic 
length [S/(mw)] 1 '' 2 is taken, the strongly-interacting 
BEC shoots out rapidly in all directions, thus producing a 
shock wave, depletion, and rarefaction waves around the 
cylindrical symmetry axis. However, the total density, 



integrated in the directions perpendicular to the symme- 
try axis, always retains a bell shape, even for very large 
interaction. 

In conclusion, we have shown that the free expansion of 
a Bosc-Einstein condensate reveals novel and interesting 
nonlinear effects, not achievable with classical viscous flu- 
ids, and which are awaiting experimental investigations. 

The authors thank Paolo Di Trapani, Luciano Reatto 
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